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Cyclic (rock-paper-scissors-type) population models serve to mimic complex species interactions. 
Focusing on a paradigmatic three-species model with mutations in one dimension, we observe an 
interplay between equilibrium and non-equilibrium processes in the stationary state. We exploit 
these insights to obtain asymptotically exact descriptions of the emerging reactive steady state 
in the regimes of high and low mutation rates. The results are compared to stochastic lattice 
simulations. Our methods and findings are potentially relevant for the spatio-temporal evolution of 
other non-equilibrium stochastic processes. 
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Stochastic interacting particle systems are a fruitful 
testing ground for understanding generic principles in 
non-equilibrium dynamics. Unfortunately, the treatment 
of such processes is marred by the absence of detailed 
balance so that the insight one has gained by analyti- 
cal means is not yet satisfactory and only few systems 
have been solved exactly [TJ [2] . Some of them serve as 
a paradigm for very complex biological and sociological 
systems. An example is the contact process, which de- 
scribes the outbreak of an epidemic, displaying a phase 
transition from an absorbing to an active state as the 
rate of infection is increased [3] . Another famous process 
is the voter model, caricaturing opinion making. It is 
proven rigorously that on a regular lattice there is a sta- 
tionary state where the two "opinions" coexist, so long 
as the dimension is larger than 2, such that the random 
walk is not recurrent [H H] . Extensive studies have also 
been conducted on the coarsening dynamics of coalescing 
or annihilating particles, both for diffusional motion and 
ballistic motion of the particles [SHE]. In this context, 
much work was devoted to the long time behavior of the 
average domain size, which as a function of time typically 
displays scaling. 

Frachebourg et. al. [3 [7J have studied the coarsening 
dynamics of a model known as the Rock-Paper-Scissors 
game (RPS), one of the most widely studied prototype 
models for biodiversity [9lUlj. displaying cyclic domi- 
nance between its three agents. In this Letter we study 
the influence of mutations on this model. An integral 
part of evolution, mutations have been posited to pro- 
mote biodiversity in microbial communities |12) . We will 
argue that the RPS is a natural framework for a nonequi- 
librium version of the Ising-Glauber model, which at zero 
temperature amounts to an annihilating random walk. 
While previous studies have addressed coarsening and the 
transition to an absorbing state, we focus on the descrip- 
tion of the stationary reactive state at finite "tempera- 
ture", i.e. interfaces between domains are created at fi- 
nite mutation rate. In the Ising-Glauber model the inter- 
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Figure 1: Illustration of the one-dimensional Rock-Paper- 
Scissors game with mutations and the passing to its dual 
description. The configuration of the lattice is given at sub- 
sequent points in time, resulting in a two-dimensional space- 
time plot. A mutation C — > B occurs somewhere between 
t — 2 and t — 3. The dual picture is characterized by inter- 
faces moving left (L) or right (R). 

faces perform a random walk, whereas for the RPS they 
drift left or right and even move ballistically in a certain 
regime. Since the coarsening dynamics is counteracted 
by the creation of interfaces, the system evolves into a 
non-trivial stationary state. For very large and very low 
mutation rates, equilibrium turns out to be only slightly 
broken. Discriminating between two types of mutations, 
we can thus obtain asymptotically exact descriptions for 
the average size of the domains in the stationary state. 
As the final arbiter of the validity of our arguments we 
employ stochastic lattice simulations. 

On a one-dimensional integer lattice {1, . . . , S} of size 
S, the RPS can be defined by the following cyclic- 
dominance reaction equations for nearest neighbors. 

AB^AA, BC ^BB, CA^CC, (1) 

i.e. paper (A) covers rock (B), rock crushes scissors (C) 
and scissors cuts paper. Here we presuppose left-right 
symmetry such that, for instance, A can invade a neigh- 
boring B to its left or right and we consider a Markov 
process in continuous time with sequential updating. Un- 
less otherwise stated, we look at the symmetric case 
ta = fB = rc and set — 1 to define the timescale. 
These equations have been studied in detail in [5J [7] . In 
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particular it was shown that, starting from some random 
distribution, the species organize in domains that un- 
dergo coarsening until finally — providing the system size 
is finite — one species takes over the whole lattice. 

In addition to the above reaction scheme for cyclic 
dominance we allow for mutations, 

A^B^C ^A, A^C ^B^A, (2) 

where we discriminate mutation cycles to the respective 
"prey" with rate fi r and to the respective "predator" with 
rate fj,i, both of which evidently conserve the cyclic sym- 
metry. The mutations counteract the coarsening of the 
system and ensure a reactive steady state. However, for 
low mutation rates, which we shall focus upon, one still 
expects the system to organize in large clusters separated 
by interfaces. Thus it is adequate to utilize the so-called 
dual description, obtained by representing the interfaces, 
i.e. the walls between the domains of one species, by 
particles (denoted L or R) and two consecutive spots 
occupied by the same species by empty sites 0. This 
mapping is illustrated in Figure [T] There are left- and 
right-moving interfaces, R and L respectively. Without 
mutations their number is bound to decrease when they 
interact: in LL — > R0 (and analogous for left and right 
interchanged, RR — ► 0L) two interfaces of the same kind 
turn into one that moves in the converse direction and in 
RL — ► 00 one has pair annihilation. 

As a starting point we derive the mean-field rate equa- 
tions. Let P(R), P(L), P{0) be the probabilities of find- 
ing an R, L or 0, respectively, at some site i. One 
then assumes the system to be well mixed, i.e. the occu- 
pancy of sites to be uncorrelated. Defining /i := /i; + fi r , 
one obtains, P(R) = -2P(R) 2 - 2P(R)P(L) + P(L) 2 + 
fj,(P(0)+P(L)-2P(R)). Solving for the stationary 
state, where P(R) = P(L), the interface density becomes 

n := P{R) + P{L) = ^(4/3)// + — £t , (3) 

increasing sharply oc ^/Jl for small rates and saturating 
to 2/3 for large rates. Notice that \x T and /j/ are treated 
on equal footing in contrast to our results below. 

Before we proceed with the RPS, we point out an 
analogy to the Ising model. It can be verified that the 
two-particle version of our process (i.e. AB AA, 
BA — ^ BB and A — > B — > A; pi and fi r obviously cannot 
be distinguished here) has been proposed by Glauber as a 
way to study the dynamic effects of the one-dimensional 
Ising model, with, say, A corresponding to "spin up" 
and B corresponding to "spin down" [T31 Q3]. After 
expressing the energy in terms of the nearest neighbor 
sum E({s}) = —JJ2 s k s h where Sk is 1 for "spin up" 
and —1 for "spin down" and J is a coupling constant, 
the temperature T is related to the mutation rate /j by 
+ fi) = 1 — tanh(2fcJ/T). fj, is small in the low 
temperature regime and large in the high temperature 
regime. Thus, we may think of the mutation rates \xi 



and \x T for the RPS as temperature-like parameters. At 
fast mutation rates, or high temperature, the system be- 
comes rather uncorrelated, and therefore mean- field ^ 
makes for a good approximation. 

Let us now try to comprehend the RPS for very low 
mutation rates /J = fj,i + ji r . Comparison with the one- 
dimensional Ising model suggests that at n = 0, corre- 
sponding to zero temperature, the system displays a crit- 
ical behavior with the correlation length going to infinity. 
In the following we show that this is indeed corroborated 
by scaling arguments as well as stochastic simulations. 

First, let us restrict ourselves to the regime \i r = and 
/ij <C 1. The interface density is low and therefore the 
single most probable mutation occurs on two adjacent, 
vacant sites 00 — ^ LR. In the particle picture this is 
achieved by, e.g. , AAA AC A. Notice that the mu- 
tation induces a predator in a — typically large — domain 
of prey, where it can spread subsequently. Hence the in- 
cidence has strong impact on the system. In the dual 
picture this is expressed in the fact that the pair LR, 
unlike RL, can separate, e.g. 

0000 0LR0 -4 L0R0 -i> L00R —>■.... (4) 

Consider what happens next to the, say, R interface. It 
moves to the right from site to site with rate 1, until it 
meets and reacts with some other interface, e.g., 

R00L 0R0L 0RL0 0000 . (5) 

It is crucial to note that diffusion becomes negligible 
when the particles are far apart, since their directional 
motion is described by a Poisson process, whose mean 
square displacement a{t) = \fi becomes small relative to 
the average distance (x(t)) = t it has traveled. Therefore, 
in our regime one should think of the particles as moving 
ballistically. For our scaling argument, we partition the 
lattice in cells of size b and consider the dynamics from 
this coarse-grained point of view. Empty cells become the 
new vacancies 0, cells that contain exactly one interface 
the new R or L, respectively. Since the lattice is supposed 
to be sparsely populated, we disregard the unlikely case 
of cells containing more than one interface. A mutation 
pi now occurs with 6-fold rate, since the whole cell is at 
disposal. We rescale time by a factor of b, so that the 
velocity of the (ballistic) interfaces is unchanged. This 
implies a rescaled rate /ij = 6 2 /_t; — one factor b for rescal- 
ing space and another one for rescaling time. The den- 
sity evidently becomes 6-fold, n(/ij) = bn(pi) and thus, 

1/2 

ri = £/ fj, t for an infinite lattice. This result is indeed 
validated by our numerical simulations (Figure [2]). 

At this point, we remark that for the symmetric case 
ta, tbi fc — 1 the stationary state can be solved exactly, 
leading to stf = V2 [15]. Here we only illustrate the 
underlying physics by the following heuristic argument. 
For symmetric rates, reactions of the type RR — > 0L 
are negligible, because it takes an interface much longer 
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Figure 2: When r^, rg, rc = 1, the rescaled interface density 
n(fj,i,fi r — Q)/^/Jm approaches the law \/2 — (3/4)2 3 ^ 4 fi 1 / 4 
(red line) as fii becomes small. For large mutation rates the 
data is described well by the mean-field curve ([3]) (blue line). 
Notice that for \xi > 1 we have plotted the density n without 
rescaling. 

to catch up with an interface of the same kind than to 
crash into a different kind of interface, which it can meet 
halfway, as it were. Hence P(RR) — 0, where P(RR) 
stands for the probability of finding two interfaces next 
to each other. We suppose that otherwise the system is 
uncorrelated, in particular P(RL) = P(R)P(L). Then, 
up to terms of the order \i\ and due to the symmetry 
between R and L one has the master equation, 

P{R) « W P(00) - 2P{RL) a w - 2 [P{R)f . (6) 

Solving for the stationary value yields stf = \/2. 

To motivate the first correction to this result, notice 
that two interfaces of the same kind move diffusively rel- 
ative to each other, with diffusion constant 1. The diffu- 
sional lengthscale associated to the average survival time 
1 /n is <l/n (for small mutation rates) . The prob- 

ability P n (t) that the pair of i?-intcrfaccs is intact after a 
time t will just be the probability that they have not yet 

interacted diffusively, ^= t ds e~ s2 (see, e.g., [T^]). 

times the probability that the right R has not yet crashed 
into an L, which is given by an exponential distribution 
exp (— nt). In the last equation we assume that the sys- 
tem is uncorrelated, so that in every time step the right 
R crashes into an L with probability n. The probability 
that an R is created at a distance x to another R is n/2 
when x <C 1/n. Upon integrating over x and t one finds 
that the probability of an R interacting diffusively with 
an R on its right is y/n/2. If we simply subtract these 
"failed attempts" of creating it from the mutation rate 
fil — > — vW2), we obtain — n/2) • 1/n = n/2 or 
n « 727^(1 - s/ri/A) « - (2 3 / 4 /4) /if /4 . More rig- 

orous analysis, taking into account the difference in the 
time of survival of the left i?s, with and without diffusion, 
yields an even larger correction 

n^^-h^^, (7) 



Figure 3: For small /j, r , the rescaled density converges to 
n(ni = 0, /i r )//i r — 2 + 3/2^/zv (red line), when rA,rB,rc = 
1. The mean- field result (blue line) and the generalized mean- 
field result for clusters containing two sites (turquoise line) is 
also shown. 

(see supplementary material) . Our result is in agreement 
with the numerics (Figure [2| . 

Let us proceed to discuss the case m = and /i r « 1. 
Since the system is coarse grained, the major part of the 
mutations will result in one prey in the middle of large 
domains of predators. For instance, AAA — ^ ABA. Ev- 
idently, this configuration is rather unstable and one ex- 
pects that in most cases the cyclic dominance reactions 
reestablish the original state, that is B turns to A. In 
the dual picture, this translates into the creation of a 
pair RL, which in most cases annihilates quickly, 

00 -^RL ^>00 . (8) 

Owing to the longevity of its products, one also needs to 
take into account that a second mutation may occur, 

00-^RL-^LR, (9) 

effectively leading to 00 >LR. Just as above a pair 

LR is produced, but this time the reaction is mediated 
by two [i r mutations instead of one fii mutation. In the 
particle picture this means that the prey B in a domain 
of A may be turned into the predator C by a second mu- 
tation. The former ^ implies a contribution of fj, r to 
the interface density. The latter ^ leads to the same 
dynamics as in the fii case and to another term of mag- 
nitude \/2 (/xp/2) = fi r , i.e. to lowest order n = 2fi r . 

For the leading correction, in addition to reactions of 
the type RR — > 0L, one needs to treat instances of mu- 
tations when there exactly one interface around, say 

R0-^LL, (10) 

e.g., ABB ACB. Similar reactions can occur, when 
there is a mutation nearby an interface, for instance, 

i?000 R0RL ± 0RRL ± 00LL . (11) 
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An analysis analogous to the pure //; case yields an overall 
positive contribution |15j . 



1\x r 
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Figure [3] confirms this behavior. Again mean-field is an 
excellent approximation for large mutation rates. As 
mutations become less frequent, equation ([3| provides a 
gross over-estimate of the interface density, because the 
approach cannot keep track of the large amount of pairs 
of RL that annihilate quickly. This can be amended by 
a generalized mean- field approach |17) , where the master 
equation for clusters of N adjacent sites is considered. A 
truncation in the hierarchy of probability distributions 
yields a closed set of differential equations, which can 
be solved numerically. For clusters of size 2 one already 
retrieves the right scaling law n ~ fi r (/i r <C 1). 

Again, we explain how a scaling analysis helps us un- 
derstand the behavior of the density n{ji r , /i/ =0). This 
will also extend our results to more general processes. We 
partition the lattice in cells of size b. Then the probabil- 
ity for a cell to contain a pair RL, which are created and 
destroyed according to reaction ([8]), becomes b[i r while 
the rate to create a pair LR out of RL ^ evidently re- 
mains [i r . We rescale time by a factor b so that the veloc- 
ity of R and L, measured in the average number of cells 
they traverse in unit time, stays one. Now the right-hand 
side of reaction ([9| occurs at a rate bfj, r , while the prob- 
ability of finding a pair LR remains bp, r . This implies 
fi' r = bjjL r , n(/jf r ) = bn, whereby we conclude n = SS[i r . 

Up to now we have considered symmetric rates, where 
reactions governed by diffusion give rise to the correction 
terms in ( 7|12 ) which do not fit in the scaling. Now con- 
sider what happens if ta^b^c are not identical, for in- 
stance if ta > tb — rc- Suppose the pair RR stands for 
CAB. In this case the two Rs no longer move diffusively 
relative to each other but rather the right R drifts away 
and can escape the left one, resembling a ballistic motion 
of the right R away from the left one with average relative 
velocity ta — rs- Here we can apply the above scaling 
argument and we expect a contribution to the interface 
density that obeys the scaling n ~ [i r as confirmed in 
Figure |4j This plot also indicates that our findings can 
be generalized to [ii 7^ and fi r ^= 0, where in analogy to 
our above arguments one derives n(fii, fi r ) = y^^^J 
for some scaling function <f>. We remark that exact cal- 
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Figure [4] shows that even 



culations yield 0(0) — 
for ta = 5 this is a good approximation 

In conclusion, competition between coarsening dynam- 
ics and mutations in our model leads to a reactive station- 
ary state characterized by an interplay between equilib- 
rium and non-equilibrium processes — indeed one can pin- 
point exactly which reactions take the system away from 
equilibrium. It was crucial to discriminate two types of 
mutations, /i; and /i r , the effect of the latter being negli- 
gible when the two rates are comparable and small. Both 
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Figure 4: Collapse of the data onto a scaling function. It 
is demonstrated that our arguments hold true also for asym- 
metric rates, here ta — 5, rs — rc — 1. The data points are 
Mi = 2- 10 (D), 2" n ( ), 2" 12 (A), 2" 13 (0). 



for the high and for the low mutation rate regime we have 
retrieved asymptotically exact results which we expect to 
be quite robust to a wide rage of variations, e.g. relaxing 
the constraints of perfect symmetry. 
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CALCULATION OF THE LEADING CORRECTION TO THE DENSITY 



In this supplement we show how to calculate the leading correction to the law n = y/2m for the interface density, 
when m <C 1 and u r = 0. Let n denote the density of interfaces that will ultimately annihilate ballistically, i.e. by pair 
annihilation RL — > 00. We expect that n ss \J1\n. In a first step, we calculate the ratio of, say, Ra that annihilate 
diffusively, i.e. via the reaction RR — > 0L. Suppose an R is created x sites to the left of the nearest further R to the 
right (there may only be empty sites and Ls in between them). For instance, if there arc no Ls in between, which is 
the most important case, this looks like 

R0 R. 
x sites 

Relative to each other the two interfaces move diffusively, with diffusion constant 1, and thus the two Rs might 
crash into each other through this diffusional motion. The expected lifetime is of the order 1/n, so that the relevant 
diffusional distances are of the order 1/y/n. Therefore, we may assume x -C 1/n in the following and neglect the 
possibility of an L in between the two interfaces. 

The probability P n (t) that the pair of i?-interfaces is intact after a time t equals the probability that they have not 

yet interacted diffusively, P<i(x,t) = ^= f^ 2 ^ds e~ s2 , times the probability that the right R has not yet crashed into 
an L, which is equal to Pb{t) = exp (— nt). In the last equation we assume that the system is uncorrelated, so that 
in every time step the right R crashes into an L with probability n, giving rise to an exponential distribution with 
parameter n. Furthermore, we note that when x <C 1/n the probability that an interface is created at x is n/2 w n/2 
(since the density of R is n/2). Thus the probability that a particle annihilates diffusively becomes 
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where 1/ \/n <C c <C 1/n. For uncorrelated interfaces, the average time until a particle annihilates ballistically is 1/n. 
Multiplying the input rate of the interfaces, 2pi, by the probability that there is no diffusive interaction, 1 — \fn/2, 

the density of particles that annihilate ballistically becomes n = 2pi (l - Vn/2) 
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n is not the interface density that we are looking for. To find n we need to add the contribution of the particles that 
annihilate diffusively. To do this we look at the average time of survival t(x) of an interface that is created at x. We 
also introduce the quantity Tb(x), which denotes the average time of survival for truly ballistic interfaces. We know 
that if diffusion can be neglected, the average time of survival is 1/n, i.e. 



(n) 



Jo 



dx p{x)rb{x) = 



n 



where p(x) denotes the probability that an R is created at x. 

Without diffusion, the average time of survival of an R that is created at x <C 1/n would be 2 • 1/n since the right 
R must be annihilated before the left one can be destroyed. Together with diffusion one has instead 



t(x) = 



dtt 



(Notice that upon ballistic annihilation of the right R at time t, the left one lives on for 1/n units of time on average, 
whence the factor t + 1/n.) For x » 1/y/n diffusion is negligible, t(x) = Tb(x), justifying the boundary c in the 
following integral. Again we remark that when i<1/b then p(x) — n/2 rts n/2. Thus 



where l/\/n «c«l/n. Therefore, 

M-i- 

n 

Expressing n in terms of m we find for the interface density 

n = 2^(r) w 



